function grad_jt = f_calc_gradient( ...
    theta,      ...
    Pi,         ...
    x_jt,       ...
    xnonlin_jt, ...
    dfull_jt,    ...
    group_ids,   ...
    share       ...
)

% Initialize the result matrix
grad_jt = zeros(size(group_ids, 1), length(theta) + length(Pi));

% Gradient calculation for the fixed effects
for i = 1:size(grad_jt,2)
    if i <= length(theta)
        product = share .* x_jt(:, i);
    else
        index = i - length(theta);
        product = share .* xnonlin_jt(:, index);
    end

    product_wid = [group_ids product];
    [Group, ~] = findgroups(product_wid(:, 1));
    sum_by_group      = splitapply(@sum, product_wid(:, 2:end), Group);
    expanded_sum_by_group = sum_by_group(Group, :);

    if i <= length(theta)
        dtheta = share.*(x_jt(:, i) - expanded_sum_by_group);
    else
        index = i - length(theta);
        dtheta = (dfull_jt.*share).*(xnonlin_jt(:, index) - expanded_sum_by_group); 
    end

    grad_jt(:,i)    = mean(dtheta, 2);
end

end
